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Abstract 

We introduce a proposal to modify Einstein's equations by embed- 
ding them in a larger symmetric hyperbolic system. The additional 
dynamical variables of the modified system are essentially first inte- 
grals of the original constraints. The extended system of equations 
reproduces the usual dynamics on the constraint surface of general 
relativity, and therefore naturally includes the solutions to Einstein 
gravity. The main feature of this extended system is that, at least for 
a linearized version of it, the constraint surface is an attractor of the 
time evolution. This feature suggests that this system may be a useful 
alternative to Einstein's equations when obtaining numerical solutions 
to full, non-linear gravity. 

1 Introduction 

Over the past decade, computer power has increased to the point that simu- 
lations of two- and even three-dimensional general relativity are now feasible. 
These simulations, which assume little or no symmetry of their generic field 
configurations, at first seemed to represent straightforward generalizations 
of simpler one-dimensional calculations. However, attempts to perform the 
higher dimensional simulations have revealed a variety of unexpected fea- 
tures which limit accurate simulations to a rather short time interval. One 
such feature, which is believed to be a major source of numerical error, is 
that numerical time evolution generates a rapidly growing violation of the 
constraint equations. In this paper, we propose a system of dynamical equa- 
tions wherein the evolution naturally remains close to the constraint surface. 
Although the most obvious application of this approach is to numerical sim- 
ulations, it may prove useful in other branches of general relativity as well. 

As is well known analytically, the time evolution predicted by the exact 
Einstein equations is such that the constraint equations are satisfied on each 
time slice when they are satisfied by the initial data. Geometrically, the evo- 
lution vector field is tangential to the constraint submanifold, implying that 
solutions to the complete set of equations are insensitive to properties of the 
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evolution field in the vicinity of the constraint surface. In discrete approxi- 
mations, on the other hand, the notion of tangency is approximate, as is that 
of the constraint surface itself. As a consequence, the numerical evolution 
becomes sensitive to possible instabilities of the constraint submanifold and 
numerical solutions are, in general, carried away from it exponentially fast 
with time. Even in case one were able to construct a code whose discretized 
vector field were exactly tangent to a discretised version of the constraint 
submanifold, the same problem would be likely to arise, as numerical errors 
on the initial data would prevent a start of the time integration exactly on 
the constraint submanifold. 

As demonstrated in ||Cho9UI, evolution schemes can be constructed in 



such a way that the violation of the constraints has the same convergence or- 
der as the scheme itself. This property, which is in the meantime a standard 
requirement for evolution schemes, implies that the choice of an appropri- 
ately fine grid is sufficient to satisfy the constraint equations at any given 
time with arbitrary accuracy. However, since the violation of the constraints 
grows very quickly with time, the utility of grid refinements to reduce con- 
straint violations is very limited, especially in two- and three-dimensional 
calculations. 

In the so-called constrained evolution schemes one, attempts to solve this 
problem by isolating two sets of variables in Einstein's equations. One uses 
evolution equations to evolve one set and determines the variables of the 
other set by solving constraint equations on each time slice. This method 
has frequently been used in one-dimensional simulations where, on the one 
hand, it is easy to split the variables into dynamical and longitudinal ones, 
and where, on the other hand, the constraint equations are ordinary differ- 
ential equations along a space-like direction. However, in two- and three- 
dimensional simulations with space-like hypersurfaces as time surfaces, the 
elliptic character of the constraint equations makes it expensive in computer 
time to solve the constraint equations on each time slice. [] Furthermore, this 
approach does not guarantee that the complete set of Einstein's equations 
is solved. Since only a subset of the variables is determined by evolution 
equations, some of these equations remain unused. The problem is, there- 
fore, shifted to the preservation of the unused evolution equations, which, as 
shown in | Det87 |, is a problem of similar nature. 

In order to guarantee a good approximation to the complete set of field 



1 Not to mention the problems arising in the treatment of the grid boundaries. 
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equations, it is, therefore, necessary to analyze the behavior of the evolution 
vector field in a whole neighborhood of the constraint submanifold. Away 
from the constraint submanifold the evolution field is not uniquely deter- 
mined as field configurations violating the constraint equations are physically 
not relevant. Hence, the evolution vector field can be modified in an arbitrary 
way, as long as its values on the constraint submanifold remain unchanged, 
and as long as the modified field continues to be strongly hyperbolic, so that 
the Cauchy problem is well posed in a whole neighborhood of the constraint 
surface. 

Of particular interest are modified equations for which the constraint 
submanifold is asymptotically stable, because for equations with this fea- 
ture, sufficiently accurate codes are expected to generate solutions which 
remain close to the constraint surface, and which, therefore, would represent 
improved approximations to Einstein's equations. 

Modifications of the evolution vector field have previously been studied. 
However, in general these preserve the time reversal symmetry of Einstein's 
equations, which implies that modifications of this type cannot have the de- 
sired properties. Time reversal symmetry implies that if the evolution field is 
such that a solution to some initial data in a neighborhood of the constraint 
submanifold approaches the constraint submanifold during time evolution, 
then the solution to the time reversed initial data will asymptotically be car- 
ried away from this submanifold. Thus, without a modification of Einstein's 
equations which breaks the time reversal symmetry, the best one can expect 
to achieve is a set of equations for which the constraint submanifold is sta- 
ble, but not asymptotically stable. However, stability of the undiscretized 
equations is not sufficient for numerical simulations, since spurious solutions 
to the discretized equations can grow very rapidly even for stable systems. 
To eliminate the impact of such solutions, it is, therefore, necessary that the 
constraint manifold is an attractor for the time evolution. 

The above-mentioned modifications of Einstein's equations are, in gen- 
eral, obtained by including dynamical quantities which are proportional to 
the constraint expressions. An alternative argument showing that exten- 
sions of this type cannot lead to an asymptotically stable constraint surface 
is the following: Since the constraint equations are of the same order as 
the evolution equations, their inclusion affects mainly the principal part of 
the evolution equations, whence the freedom remaining after requiring that 
these terms do not destroy strong hyperbolicity is very limited. Thus, such 
extensions ensure the well posedness of the problem, but not the asymptotic 
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stability. This can only be obtained either via modifications of the lower or- 
der terms or the addition of higher (second) order terms, that is, by including 
damping or diffusion terms. 

In the next section, we propose a modification of Einstein's equations 
which includes new dynamical terms proportional to certain first integrals of 
the constraint expressions, rather than to the constraints themselves. The 
dissipation, that is the time asymmetry is not of the diffusive type, 2 and is 
built into the definition of these integrals. 

We show that the Cauchy problem for the resulting new system, which 
we call the A-system,Q is locally well posed. We also prove that if the 
constraints are initially satisfied, and if their first integrals initially vanish, 
then the A-system provides solutions to Einstein's equations. Moreover, for 
initial data sets for the A-system, which are sufficiently close to the constraint 
submanifold and sufficiently close to zero, respectively, we suspect that the 
solutions asymptotically tend to solutions to Einstein's equations. 

In section [| we give support to our expectation by proving that the 
linearized extended system is asymptotically stable, thus showing that in the 
linearized case, the constraint submanifold is indeed an attractor for the A- 
system. In section f|, we discuss further expectations in connection with our 
proposal. 



2 The A— system 

In this section we spell out our proposal for a modification of Einstein's equa- 
tions with an asymptotically stable constraint submanifold. For definiteness, 
we choose the symmetric hyperbolic system introduced by Frittelli and Reula 
in |[FR94|| , which corresponds to the parameters a — /3 — — 1 in ||FR96|| . Al- 



though the full equations (with the non-principal part terms added) are given 



in [(5te98|| , we repeat them for completeness. 

In the version of Einstein's vacuum equations chosen, the system is given 
by the following set of dynamical equations (where Latin indices run from 1 
to 3): 

■ h H = N n h lJ , n +QVh (2P ij - Ph l] ) - 2h n{i N j \ n , (1) 

2 Onc could also introduce diffusive dissipation, but this would significantly reduce the 
allowed maximal time step in explicit discretisation schemes. 

3 The name is a remnant of the way the system was originally guessed by Brodbeck and 
Hiibner, namely by a formal application of Lagrangian multiplier techniques. 
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W\ = N n M^ k ,n + QVh(P l] ,k-2h {l P j)n rn) 

+ QVh(^P l] M k - PM l \ + Q^P^Q, 



- 25 k ^\h^h mr h ns P mn M rs q - 2M^ n P mn h pm 



+ ^P j)n M n - hi )n PM n 

+ 0N n mk -h^Ni\ nk -2N«, n M^ k + N m )k M^ m , (2) 
pij = N npij m +Q ^h (h mn M^ m , n - 2/i n < i M'> fc>n ) 

+ QVh(^h np h m{i M^ n k M kp m - K k h? n h rp h sq M r \M™ n 

+ hj k h jn M k M n + 2M nk k M ij n - 2M ik n M jn k 

- 2h mn h kp M im k M jn p - 2h n{i M j)k k M n + M ij n h nk M k 

- Q- 1 [h tk V n Q, kn +2M k (\V )n Q, k -M^ m h mk Q, k 

- (h kn Q, kn +2M km m Q lk ) ] + 2P lk h kn P^ - ^-PP^ 



+ V j (±P 2 -h mr h ru F mn F r ')^-2P k ( i W\ k . 



(3) 



Here, h l i is the inverse intrinsic metric of the spacelike hypersurfaces S t , 
pij ._ fcij _ foijfc d erio tes a linear combination of the extrinsic curvature k l i 
of the slice and its trace k, and M tj k := \ (h lj jk — h %] h rs h rs :k ) represents a 
linear combination of spatial derivatives of the inverse intrinsic metric. The 
functions P and M k are abbreviations for h^P^ and hijM tJ k , respectively, 
and Q and N l are arbitrary given functions fixing the gauge degrees of free- 
dom. 

This evolution system is symmetric hyperbolic with respect to the inner 
product 

(hlPl\M? k \H e \hlP?,M? k ) : = 

»Mi j k M? n i} dE, (4) 

where denotes an Euclidean flat metric on the hypersurface E t . It is 
completed by the following set of constraints equations: 

C = 0, C l = 0, C^ = 0, (5) 
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where 



C 




- h mn h rs h™M mr p M n % - hx mn h rs P mr P ns + -^P 2 , (6) 
P lk , k -2h mn M m k P nk - l -h lk PM k 




The first two constraints are the scalar and the vector constraint of Ein- 
stein's equations, that is the time-time and time-space components of the 
Einstein tensor for a given 3+1 decomposition of space-time. The third is 
the statement that the tensor M %3 k is a linear combination of spatial deriva- 
tives of the 3-metric. 

To solve the initial value problem of general relativity in this approach, 
one prescribes an initial data set (Hq , Pq , Mq k ) at t — which satisfies the 
constraints equations and subsequently solves the above evolution equations. 
Symmetric hyperbolicity of the evolution system implies that a unique local 
solution does exist. 

By taking a time derivative of equations and using ([[]-§]) to elimi- 

nate time derivatives in favor of spatial derivatives, the following evolution 
equations for the constraints are obtained: 



where ". . ." represent undifferentiated terms which are linear in the constraint 
quantities and at least linear in the variables P % i and M lJ k . 



C 



N n C, n +3QVhC\ k + ... , 

N n C\ n +QVh (h lk C )k +h rs C k [rMs + h ls h kl h mn C mn [sMl ) 



(9) 



+ ... , 

N n C ij k , n - 2QVh (25 k {i C j) - h ij h kl C l ) + . . . , 



(10) 
(11) 




4 Onc could also consider the constraint C l j := C lk j k , which still makes the constraint 
system symmetric hyperbolic and produces a smaller number of extra fields. 
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By taking a time derivative of (|T2D, we obtain 

C« fcJ = iV™^ H , n - 2Qv^ (& (< C J '>,, S^C j \ k )+..., (13) 

and by plugging (|i~2"D into (|10|), we see that the evolution equation for C 1 can 
be rewritten as 

C = N n C\ n +QVh (h' lk C, k +h rs C k rk , s ) + ... . (14) 

The constraint quantities C, C\ C u k , and C^'jy thus propagate according 
to the first-order system of equations consisting of (|), (|TTD, and (|13Hl4"|) , 
which is symmetric hyperbolic with respect to the following inner product: 

^1) ^1 k, W kl | DC | ^2, (-"25 ( - y 2 fc, <-^2 fci/ - _ 
J St ^ O 

V'C^A™^}^- (15) 

Uniqueness of solutions to this system implies that if the constraints are ini- 
tially satisfied, then the exact evolution equations preserve them. When, as 
in numerical simulations, the constraint variables initially are not precisely 
zero, then the corresponding solution is, in general, carried away from the 
constraint surface during time evolution. However, since the evolution equa- 
tions for the constraint variables are symmetric hyperbolic, the violation of 
the constraints becomes smaller when the constraints initially are satisfied 
with better accuracy. 

In order to obtain a system with an asymptotically stable constraint sub- 
manifold, we propose a modification of Einstein's equations, which is inspired 
by the behavior of dissipative systems, where a transient eventually is dis- 
sipated away as the system settles down. We extend the set of dynamical 
variables by considering the following "time integrals" of the constraint vari- 
ables: 



A = a C-A)A, (16) 

X = ai C - /3 1 A i , (17) 

X ij k = a 3 C ij k -p 3 \ ij k , (18) 

\ ij M = ot^ja-fa&u, (19) 
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where the tensor-valued A-variables are assumed to have the same symmetries 
as the corresponding C-variables, and where «j 7^ and $ > are constants. 

The equations (pl6|-|T9|) represent evolution equations for the A-variables 
which in terms of the fundamental variables (h^, P y , M'^) are given by 

A = a (-M fc \, fc + h pq MV k n M« n k - M k \M k + l -h kn M k M r 

- /?oA , (20) 
A* = a l ( K P t \ k -2h mn M im k P nk -h* k PM k ^-[3 l \\ (21) 

X'\ = a 3 (2M i \-h ij h rs M rs k -h ij , k )-P 3 X ij k , (22) 
X ij kl = a i (2M l \ kA + 2M i \ k M l] )-f3 i Xhi- (23) 

In the present form, the combined system ([T| |3|j20| ^3f) is not symmetric hy- 
perbolic, since the equations (p0|-p3|) involve spatial derivatives of the vari- 
ables (0 , P % \ M y jt), whereas the equations (0-0) do not contain A-variables 
at all. However, by adding terms containing first derivatives of the A-variables 
it is possible to bring the system (|T|-|3,20-F23D into a symmetric hyperbolic 
form, 

& = a 3 h mn X ij myn + N n h ij m + source terms , (24) 
W\ = 2a i h lm X ij kl>m -a Q S k ^ l X, l 

+ N n M ij k>n + QVh (pv, k -25 k ^P^ n , n ) 
+ source terms , (25) 
p« = a 2 /i^A^ 

+ N n P l \ n +Qv^ {h mn M^ m>n - 2h n ^M^ k k>n ) 

+ source terms . (26) 

By construction, the "A-system" (|20|-p6|) is symmetric hyperbolic with re- 
spect to the inner product 

{hi ') P\ ■> Mi k , Ai, A^, Xi fc, Xi k i I I h 2 , P 2 , M 2k , ^2, A 2 , X 2 k , X 2k i) '■— 
[ {e im e jn h[ j h™ + e im e jn Pl j Pr + e im e jn e kl M? k M™ n l + A X A 2 

+ + CipCjqC X-f kX^r &ip&jq& & ^/fcZ^l^rs | dYl . (27) 
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The initial data for this purely dynamical set of equations consists of 
arbitrary functions 

(ho, Po, Mq\, A , Aq, Aq k, Ao ki) ■ (28) 

However, the dynamical degrees of freedom are extended by 40 A-variables. 

Clearly, for an arbitrary solution to Einstein's equations, (h%, Pg , M l ^ k ), 
the embedded field configuration (h ij , P ij , M^ k , A, A*, \ ij k , \ ij M ) := (h%, P%, 
M% k , 0, 0, 0, 0) is a solution to the A-system. Conversely, every solution to the 
A-system with vanishing A-variables is also a solution to Einstein's equations. 
Due to this property, and since the solutions to the A-system are unique, the 
A-system naturally reproduces the dynamics on the constraint submanifold 
of general relativity. 

Note that if the constraints are initially not satisfied, then, even when the 
A-variables initially vanish, the A-variables would pick up a non-zero value 
during time evolution. Hence, solutions to the A-system corresponding to 
such initial data sets would not represent solutions to the complete set of Ein- 
stein's equations. In fact, they would not even solve the evolution equations 
of general relativity. However, for constraint- and A-variables which initially 
are sufficiently close to zero, we suspect that the solutions asymptotically 
approach solutions to the Einstein equations. In the following section, we 
give analytical evidence that this conjecture could be true. 

The system is by no means uniquely "extended" , since one could still add 
non-principal (undifferentiated) terms, as long as they vanish when A = A 1 = 
A u fe = A* J \i = 0. Such terms might be useful in order to treat the strongly 
non-linear regime. Of particular interest might be to choose the coefficients 
aii and fa, which control the damping in the A-equations, to be quadratic 
functions of the basic variables (h' lj , P' lj , M' lj k ), so that the damping becomes 
stronger at points where the non-linearities intensify. 

It is fairly easy to implement similar schemes for alternative symmetric 
hyperbolic systems for the Einstein equations, as well as for symmetric hy- 
perbolic systems for other theories with constraints, like, for instance, Yang- 
Mills theories. The strategy is the same: One writes equations with damping 
for first integrals of the constraints and modifies the evolution equations such 
that the extended system becomes symmetric hyperbolic. This can always be 
achieved, because the inclusion of the new equations modifies an off diagonal 
sector of the principal symbol matrix. 
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3 Asymptotic stability of the constraint prop- 
agation 

The inclusion of the A-terms into affects, in turn, the evolution of the 

constraint quantities C,C l , C^k, and C^ki- Recalculating the time derivative 
of these, and using (p^-p6|), yields the constraint evolution equations for the 
new system, 

C = N n C, n +3QVkC\ k -2a 4 h mn X kl km>nl + 2a h mn X, mn + . . . , (29) 
C = N n C\ n +QVh(h ik C, k +h rs C ik rKs ) + a 1 h m ^ n X l \ mn + ... , (30) 

C\ = N n C ij kjn -2QVh(2S k ^C j) -h ij h kl C l ) 

+ 2a 3 h mn X lj m , nk + 2a 4 h mn (2X lj km , n — h lj h TS X rs km , n ) 
-a (25^ k h^ l X, l -h ij \, k )+ ... , (31) 

C ij kl = N n C ij kl>n -2QVh(6 k ( - i C j \ l -5 l ^\ k ) 

+ 2^4 {h™ 1 "" 1 ^ km ,nl — X V lm,nkj 

- a (S k {l h^ m X, ml W l h j)m X, mk )+.... (32) 

Again ". . ." represent undifferentiated terms that are linear in the constraint 
quantities and at least linear in (P^, M^ k ). 

The propagation of the constraints is ruled by the system of equations 
consisting of (p!6Hl9|) and ([29| |32|) , which determines whether or not the con- 
straints asymptotically "decay" to zero. The crucial feature of this system is 
that the right hand side also contains non-principal terms. Roughly speak- 
ing, the operator on the right-hand side amplifies constraint violations if the 
matrix representing its action on periodic functions has any eigenvalue with 
a positive real part. On the other hand, if all the eigenvalues have a negative 
real part, the operator induces an asymptotic decay of these violations. 

Instead of attacking the full non-linear problem as stated, which repre- 
sents a problem well beyond the scope of present analytical techniques, we 
consider the linear regime of general relativity. That is, we restrict attention 
to 3-metrics of the form = e y + with e*- 7 = 5^ and e C 1. This 
implies that the variables (P iJ , M l ^ k ) are of first-order in e, as are the con- 
straint quantities (C, C\ C^ k , C 13 k i) and the variables (A, A*, X^ k , X^ k i). Thus, 
the terms represented by ". . ." in the equations (p9]-|32D are of second order 
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in e and shall be neglected. Without loss of generality, we restrict the follow- 
ing arguments to the case where the gauge source functions Q and N l are 
constant. All arguments that follow refer to this linearized regime. 

Although we lack a proof for the non-linear case, the following considera- 
tions provide analytical evidence for the asymptotic stability of the constraint 
propagation, in particular since the full evolution equations are quasi-linear. 

For, as we believe, purely technical reasons, we adopt the following choice 
of coefficients: (3q = f3\ = /3s = f3± := (3 > and 04 = ^«o- 

Theorem 1 With the above assumptions, the constraint submanifold of the 
linearized Einstein equations is an asymptotically stable submanifold for the 
solutions to the linearized, X- extended Einstein equations. 

We partition the proof of this theorem in several lemmas: We first show 
that the initial value problem is well posed and that the solutions stay 
bounded with time. Thus, it is possible to apply Laplace transformation 
techniques, which reduce the problem to the study of the eigenfrequencies 
of the system. For these frequencies, we show that the real part is non- 
positive, only approaches zero as the wave number goes to zero, and does so 
quadratically. Then stability follows from estimates in |[KKL98|| . 

Without loss of generality, we expand the linearized dynamical fields in 
Fourier integrals of the the following form: 

X(x, t) = J X (k,t) exp(ik-x ) d 3 k , (33) 

A* (a;,*) = J A* (k,t)exp(ik-x) d 3 k , (34) 

i (35) 
C ij k {x,t) = J C il k {k,t)exp{ik-x) d 3 k , (36) 

C ij kl (x, t) = f C ij kl (k, t) exp(ik-x) d 3 k , (37) 



where k-x := kix % . 

In terms of the Fourier transformed variables, equation fl29H32"l) and (|16|- 



19) reduce to the system of ordinary differential equations given by 



A = -P X + a C , (38) 
A* = -(3X l + a x C l , (39) 
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\ ij ki = -^hi + ^-^kt, (40) 
C = ik n N n C + 3iQk m C m + V3a \ rl rm k m k l -2a \k n k n , (41) 
& = ik n N n & + iQ (k'C + k r C in rn ) - ^a^knX* + k' l k n X n ) , (42) 

bhi = ik n N n Chi - 2iQ (5 k ^h - 6fi6»k k ) 
- V3a (\ ij krk r h - X ij lr k r k k ) 

+ a \(6 k ( i k»k l -6 l Wk k ) , (43) 

and 

%ij k = _^ k + a 3 Ch, (44) 
&\ = %k n N n & k - a 3 X^ m k m k k + S l h , (45) 

where 

S ij k := - 2Q (28 k Hfi - 5 i3 C k ) - a (26« k h» - Wk k ) X 

+ V3a k m (2\v km - h ij h rs X rs km ) ■ (46) 

This system of equations naturally splits up in two subsystems, since the 
equations (|38| -4l^) couple to the equations (^4| 445|) only via the "source" term 



in (pE5|). In the following, we will first establish that the solutions to the sub- 
system (j38|- |4"2"l) ; and hence the coupling term in (fj5f), asymptotically decay 
to zero. In a second step, we consider this coupling as a given, decaying 
source, and discuss the asymptotic behaviour of solutions to the subsystem 



(EHUD- 



Lemma 1 LetTi be the space of the Fourier transformed A* J 'jy G L 2 , and let 
T> C 7i be the subspace defined by X^ ks k s = 0. Then T> is invariant under 
time evolution, and the trivial solution X J k i = is asymptotically stable for 
the evolution restricted to T>. 

Proof: Multiplying equation ( [10]) by k m , antisymmetrizing, and using that 
C l \ k ik m] = M l \ k hk m] = 0, we obtain 

X ij [k ik m] = -{3X ij [k ik m] . (47) 
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Next we note that for a function X^ki in W, the component (A 1 -^)" m ^ is 
given by (A y 'jy)H = & lj k r e rk i, where a lj = X lj [ k ik r ]S klr / (6k 2 ) . Equation (f47|) 
is, therefore, equivalent to 

C\ l h^ = -f3(X tJ ki) l] , (48) 

which proves lemma [I]. 

By direct inspection of the evolution equations, it follows that the equa- 
tion for the component of X ^i in the subspace V decouples. It is, therefore, 
sufficient to concentrate on the evolution in the space CF\ © CFq which 
comprises those functions (X, X , X ^ k i,C,C' l ,C^ki) £ L 2 for which ) \ iJ ' ki e V ± . 
Here, V 1 - denotes the L 2 complement of T> in 7i, which, as easily seen, is 
spanned by the elements X^m £ L 2 satisfying X \kik m ] = O.f] Since for the con- 
straint variable C^^i, the same property is fulfilled, C % \uk m \ = 0, this shows 
that the spaces CF X and CF C are naturally isomorphic, CF X ~ CF C = : CF. 

To simplify the notation, and to display the structure of the evolution 
equations considered more transparently, let us introduce the following op- 
erator E acting on functions v := (v , , m) in CF: 



where 



E(v) := (E(v),E l (v),E^ kl (v)) , (49) 

E(v) := V3a v rl rm k m ki - 2a vk n k n , (50) 

E\v) := - l - ai (v*k n k n + v n k l k n ) , (51) 

E ij hl (v) := -Vla,(y l hnk n h-v^ ln k n k k ) 

+a v(5 k {t k^k l -5 l il k j) k k ) . (52) 



Taking advantage of these definitions, the evolution system (|38[ - f43D restricted 
to the subspace CF © CF can be rewritten as 

, (53) 





5 For functions X ^ k i m D on ly the components along k m are non-trivial, A ^ m = 
—2knX : 'k]mk m /k 2 . This can be seen by solving X ^\kik m ] — 0, and by using the anti- 
symmetry in the lower indices of A J ^ . 
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where S and T are diagonal matrices determined by the parameters (3 and 
ccj, respectively, and where A is an operator of the form A m k m . 

In a next step, we show that the operator e is bounded with respect 
to a suitably chosen norm. To this end, we first establish the following 

Lemma 2 The operator H\ := —T~ l H c E considered as a matrix-valued field 
on the Fourier space B? is symmetric and coercive with respect to the inner 
product (u,v) := uv + CijU % v^ + ei P ej q e kr e ls u 1 ^ kiv pq TS - That is, (u,H\v) = 
(H\u, v) for all u, v G CF, and there exists a constant c > such that 
(u, H\u) > ck 2 (u, u) for all u E CF. 

Proof: We have 

(«, T- 1 ^ E(v)) - {T- l H c E(u),v) 

= -^-u (-2a vk n k n + V3a v M km k m k l ) 



+ -^u* (-^ (v 3 k n k n + fcVfc, j j r „ 



2y / 3ao 

— u (-2aoMfcnfc n + V3ct u kl km k m ki 



3ao 
1 

at i 



—v l (-^(u ] k n k n + k j u l ki) ) e i3 



+ ^e m e J „e t! 'e^^ 1 (-2V3a u mn ps k s k q + 2a u5 m p k n k q 
= 0. 

The remaining part of the proof is given in appendix [A], where we show that 
H\ is coercitive with constant c = 1/4. 

With the help of lemma 0, it is now easy to prove 



Lemma 3 The matrix-valued fields P±, 

( n F \ 

(54) 



v= 5 p ■= r 

+ ■ ' E M 



are hermitian respectively anti-hermitian with respect to the inner product 
((Ax, d), # T (A 2 , C a )> := (A 1; # A A 2 > + (d, tf c C 2 ) . (55) 
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Proof: Since S = (31, the statement for P + is trivially true. The anti- 
symmetry of P_ follows directly from lemma 0, and the symmetry of A with 
respect to H c . 

Taking advantage of lemma |3|, we now obtain the following important 
estimate for the operator P = P+ + P_: 

H T P + P ] H T = H X S + SH X = -2(3H X < -2(3^k n k n < , (56) 

where, for any Hermitian matrix M, the inequality M < means (v, Mv) < 
for all vfl 

The symmetry and coercivity of the operator H x imply that H x can be 
used to define a scalar product on a (dense) subspace T>(CF) of the Hilbert 
space CF^\ which, in turn, shows that the operator Hx = H c + H x gives rises 
to a scalar product on CF © V(CF). 

As is well known (see, for instance, |[KL89|| ), the estimate ( |56"D implies 



that for all t > 0, the operator e^ >t is bounded with respect to the norm 
defined by H?- Hence, the initial value problem for the system considered is 
well posed. Moreover, all solutions with initial data which are bounded with 
respect to this norm remain bounded for all positive times. Thus Laplace 
transformation techniques can be applied [|KL89|| , and the relevant questions 



are the sign of the real part of the eigenvalues of P, and how fast they 
approach zero as the wave number k = y/Wkl goes to zero. Hence, the proof 
is reduced to the eigenvalue problem for the operator P, 

(57) 





Then we have the following 

Lemma 4 The eigenvalues of the above system have non-positive real part 
and furthermore there exist positive constants c\ and w\ such that 

< "Ci^Tj (58) 

for all wave vectors ki 

6 Clearly, there are other possible choices of the operators S which lead to the same in- 
equality. Here we have restricted to the simplest possibility, but for practical applications, 
alternative choices might be better suited. 

7 In physical space, the relevant function space equipped with the norm corresponding 
to the above scalar product is very similar to the Sobolev space Hq . 
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Proof: From the A-rows of the eigenvalue equation, we get 



C a = (s + £)]?-%. (59) 
Using this in the C-rows, we next obtain 

(E + (s + (3)(-sI + iA)T- 1 )X s = . (60) 

Multiplying the above equation by — (T^ 1 yH c from the left and subsequently 
contracting with A s , we find the following second order equation for the 
eigenvalue s: 

(A a , H x \ s ) + (s+(3) (s(X a , (T-yHcT-'Xs) - i(X s , (r -1 )tff c .4.r~ 1 A a )) = . 

(61) 

The established properties of the involved operators imply that 

is positive for ki ^ 0, and that 

._ (A. (T-yH c AT^X s ) 

h[K)k - (\ s ,(t-^h c t-i\ s ) (63) 

is real, where ft® denotes the unit vector in the direction of k iy and k is the 
norm of ki. Thus we have for each direction of k l 

(s + (3)(s - ibk) + ck 2 = (64) 

with (3, b, c real and /3, c positive. For this equation we prove in appendix |B| 
that the real part of the roots satisfies the desired inequality, which establihes 
the result for each direction of the wave vector fcj. Using the maximal values 
of — c\ w ^ +k z on the 2-sphere of directions of k{, we obtain the final inequality. 

With this bound on the decay constants, it is now easy to prove asymp- 
totic stability for the subsystem (08|-|43|). Splitting the set of solutions in 
a part with frequencies with k < 1, and another with k > 1, the above 
bound tells us that the solutions of the higher frequency part decay faster 

c-^ t 

than e , while the decay of the solutions of the low frequency part can 
be estimated as in ||KKL98| , lemma 1 and 2 of section III]. 

We now turn attention to the second set of equations, given by (fH]) and 
(|45|), and establish the following 
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Lemma 5 Let 7i 3 be the space of the Fourier transformed (A y fc, C^k) £ 
L 2 . Then 7i 3 is invariant under time evolution, and the trivial solution 
(X lJ k,C tJ k) = is asymptotically stable for the evolution restricted to 7i 3 . 

Proof: In a first step, we discuss the equation for the component of a solution 
in the subspace 

X? 3 := { (X tJ k, C l \) e L 2 | X^ m k m = C^ m k m = } . (65) 
Taking advantage of equation (H), (|12D, and (f|4]), we obtain 

~\ is \ikn = -/3X ii lt k l] +a 3 C'\ t k ll , (66) 

which implies that the space T> 3 is invariant under time evolution. As already 
shown, the constraint variable C lJ \i asymptotically decays to zero. The dy- 
namics in T> 3 is, therefore, described by a system of ordinary differential 
equations of the form u = — u + /, where / is a given source with / — > 
as t — > oo. Since any solution to this system satisfies u — > as f - * ooJf| it 
follows that solutions in D 3 decay with time. 

It remains to discuss the complementary subspace Pg , 

Vi = { (X'h, C'h) e L 2 | X ij [k k l} = &\ k k l} = } . (68) 

For the component of a solution in this subspace, we find 

± ( A 3 \ ( -P a 3 
dt\C 3 ) " ^ -a 3 k 2 ik m N' 

where (A3, C3) := (X tJ \, C 13 k)' L G ^3", and where F %3 kk is a shorthand for the 
perpendicular component of the source term S lJ \, F ij = S lJ m k m /k 2 . Thus, 
as expected, the subspace V 3 is invariant as well. 

Since P l] and consequently C l /\k\ = P im k m /\k\ are contained in L 2 , 
equation (|46| ) implies that the same is true for F %3 \ F l i G L 2 . Furthermore, 

8 For a proof, choose T such that f(t) < e/2 for all t > T. Since the general solution 
to the above system is given by u(t) = e~*(u(0) + f e t f(i)dt), it follows that u(t) < 

e~*(u(0) + jT e*f(t)dt — ee T /2) + e/2. Hence, for a sufficiently large time to > T, the 
absolute value of the first term becomes smaller than e/2, which implies | u(t) \ < e for all 
t > t - 
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the real part of the eigenvalues of the system (j^j can, as in lemma |], be 
estimated by the inequality (|58|), albeit for different constants. Adopting a 
similar reasoning as in the previous discussion, and applying lemma 1 and 
2 of ||KKL98|| to this system, it follows that solutions in T>^ also decay with 
time. 

This completes the proof of lemma |5] and hence the proof of our main 
result. 



4 Conclusions 

In the present paper we have shown that an arbitrary system of symmetric hy- 
perbolic evolution equations with constraints admits extensions to symmetric 
hyperbolic systems which reproduce the original dynamics on the embedded 
constraint submanifold. We have given analytical evidence that the class 
of extensions proposed is sufficiently rich to contain systems for which the 
embedded constraint submanifold is an attractor of the time evolution. For 
the Einstein equations, we have constructed an extended evolution system 
for which, at least in the linearized case, this property is fulfilled. 

It is natural to expect that, by making use of techniques developed 
KKL98| , the results proven for the linearized Einstein equations can be 



in 



generalized to the regime of non-linear general relativity describing space- 
times in the vicinity of Minkowski space. However, to establish similar re- 
sults for more extended regions of the phase space of general relativity is well 
beyond the scope of present analytic techniques. 

Numerical experiences with the Navier-Stokes equations for incompress- 
ible fluids show that asymptotic stability of the constraint submanifold is 
essential for accurate results | Kre|| . For this system, techniques with a very 



similar effect have been used to include the incompressibility constraint into 
the evolution equations. On the basis of this observation, and the results es- 
tablished for linearized gravity, we suspect that the extensions of Einstein's 
equations constructed could be of interest when obtaining numerical solu- 
tions to general relativity. Numerical experiments testing aspects of this 
conjecture are in progress. 
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A Proof of coercivity 

In this appendix we show that 

(u,H x u) > \k n k n (70) 

for unitary u satisfying u lJ rs = u^\ rs \ and u % \ rs ki-\ = 0, as needed for 
lemma ^|. We treat this as the problem of extremizing the quadratic function 
of u on the left-hand side of ( [70]) under the constraint condition (it, it) = 1. 
From (|50rl52T), we obtain 



(it, H x u) + rk n k n (1 — (it, it)) = 

2 1-1 \/3 

k kfi%LVL ~\~ k k n u Tii ~\~ u ki%Lfik ~\~ u ^ rs k iiij k m uu ^ %sk kj 

3 2 2 3 

--^=uui/ s k s k j + rk n k n (l-uu- u l Ui - u 13 rs W;/ s ) , (71) 

where r is a Lagrange multiplier and where indices are raised and lowered 
with Cij. To simplify the algebra, we choose a basis in which k n = (0,0, k). 
Then u l] rs = 0, except when s = 3. Hence, 

F(u, r) := (it, H\u) + rk n k n (l - (u, it)) = 



,9/2 _ ,1 7 '_ , 1 3_ ji — 1 — iH 

k I —ww + -w m + -u w 3 + w J r3 u i:j - ~^uu i3 

--^=uu i3 l3 + r (l - wW - w*Wi - u lj rs Uij rs ^j . (72) 

The function -F(it, r) is extremized at points (it, r) where 

dF OF 
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Equation ([75]) is the requirement that u has unit length. Equation (|73|) and 
(|74l ) constitute a homogeneous linear system of equations for the real and 
imaginary parts of u. Since u cannot vanish, the determinant of the linear 
system has to vanish. Up to numerical factors, this is given by 

(2r-l) 13 (r-l) 2 (r-i) . (76) 

As easily verified, r = 1 yields the following minimal value of F(u, r) (when 
evaluated at unit u such that (|T5|)-(|T4|) are satisfied): 

F(u min , 1) = -k 2 , (77) 

from which (|7(J) follows. The other extreme values of F(u,t) are (5/3)k 2 
and (l/2)k 2 for r = 1/6, 1/2. 



B On the proof of lemma 

In this appendix we prove that the roots s± of the polynomial 

P( s ) = s 2 + s((3 - ibk) + ck 2 (78) 
are subject to the inequality 

»(«±) < -ci^T^ > ( 79 ) 

where Ci = (3c/ (b 2 + 4c) and u>i = /? 2 / (b 2 + 4c). As in the body of the text, 
it is assumed that the parameters of P are real, and that (3 and c are strictly 
positive. 

To begin with, let us rewrite the polynomial P, and the above estimate 
in terms of suitably rescaled parameter. Defining 

s = s/P, k = 2kVb 2 + 4c/(3 , b = b/Vb 2 + 4c , (80) 
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and dropping tildes, we obtain for the polynomial 

p/(3 2 = s 2 + s(l - ibk) + (1 - b 2 )k 2 /A . (81) 

The estimate for the roots in terms of the scaled parameters assumes the 
form 

*( s± > - ~TTTF' (82) 

where 7 2 := 1 - b 2 G (0, 1] . 

As easily verified, the roots of the scaled polynomial satisfy 



max{ 9fc( 2s+) , 3?( 2s_) } = -1 + | & Vl - k 2 + libk \ . (83) 
It is, therefore, sufficient to show that 



MVl-k 2 + 2ibk I < 1 - i • 

2 1 + k 2 



(84) 



To give a proof of this inequality, we first evaluate the identity 

2|Kv^| 2 = |*| +9fc(z) (85) 

for 2; := 1 - /c 2 + 2i6A; , 

2 | V^l 2 = + + (1 - A; 2 ) 

= ^(1 + k 2 ) 2 - 4 7 2 A; 2 - (1 + k 2 ) + 2 . 

Hence, 

isft^l 2 = l + (l + A; 2 )| v /l-47 2 A; 2 /(l + P) 2 - l|/2 

< 1-7^, W 



where we have used the estimate y/1 — x < 1 — x/2, which holds for x < 1. 
Taking advantage of the latter estimate once again, it follows that 

which completes the proof of our claim. 
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